$onText
DISE-2024 (Dynamic Integrated Space-Economy) model
File: DISE_2024D (Decentralized economy)
July 2025
Reference: Optimal path for orbital debris.
Bongers, A. and Torres, J. L. (2025).
$offText

*********************************************************************
*	      Simulation period
Set t	'Time periods'				/2023*2152/;
*********************************************************************


*********************************************************************
*       Time asignations
*********************************************************************
Sets
        tfirst(t)   'First period (2023)'
        tlast(t)    'Last period (T)'
    	tnotlast(t) 'All periods but last'
    	;

        tfirst(t)$(ord(t)=1)       =   yes;
        tlast(t)$(ord(t)=card(t))  =   yes;
        tnotlast(t)                =   not tlast(t);
*********************************************************************



*********************************************************************
Parameters
*********************************************************************
**	    Economic parameters
*********************************************************************
*  	    Preferences parameters
        rho	       'Intertemporal pure preference rate'		 /0.015/
        sigma	   'Relative risk-aversion'			         /1.5/
*********************************************************************

*********************************************************************    
*       Exogenous growth sources
*********************************************************************
    	gN0	    'Labor growth rate parameter'		    /0.05/
    	Nstar	'Asymptotic population'				    /10200/
    	gA0	    'TFP growth rate'				        /0.015/
    	deltaA	'TFP growth rate trend'				    /0.001/
    	gQ0	    'Satellite ISTC growth rate'			/0.030/
    	deltaQ	'Satellite ISTC growth rate trend'	    /0.005/
*********************************************************************
* 	    Initial values economic variables
	    k0	    'Initial Earth capital'				     /555.6987/         		
    	s0     	'Initial Satellite capital'			     /1.1959/  			
    	y0	    'Initial output'				         /184.65/
    	x0      'Initial value of satellites destroyed'
    	N0      'Initial population'				     /8056/
*********************************************************************
*  	    Technological parameters      
	    alpha1 	'Output-capital elasticity'			     /0.3479/
    	alpha2 	'Output-satellite elasticity'			 /0.0021/ 
    	deltak 	'Capital depreciation rate'			     /0.07/
    	deltas 	'Satellite depreciation rate'			 /0.15/ 
	    a0	    'Initial value of TFP'				      
	    q0      'Initial value of satellite ISTC'		 /1.00/
*********************************************************************
*       Launch cost
        b0      'Initial proportion of launch cost'     /0.3/
        gb0     'Initial change in launch cost'         /0/
        deltab  'Launch cost growth rate trend'         /0/
*********************************************************************
*       Physical parameters
        deltaf	'Fragments decay rate'				       /0.01/
        deltaw  'Dereclit satellites decay rate'           /0.00015/
        deltaz  'Rocket bodies decay rate'			       /0.00015/
        omega 	'Number of fragments per launch'		   /4.00/
        xi	    'Fraction of derelict satellites in orbit' /0.40/
        theta	'Collision risk'				           /1.25e-10/
        epsilonw'Fraction of dead satellite breakups'	   /0.001/
        epsilonz'Fraction of body rockets breakups'		   /0.0012/
        psi	    'Rocket bodies per launch'			       /0.60/
        phiw 	'Fragments per derelict satellite breakup' /44.60/
        phiz    'Fragments per rocket body breakup'		   /100.20/
        gammas	'Fragments per operational satellite collision'	/70/
        gammaw  'Fragments per derelict satellite collision'	/70/
        gammaz  'Fragments per rocket body collision'		    /70/
        eta	    'Number of satellites per launch'		        /13.6/
        v	    'Collision avoidance'				            /0/
        mu	    'Mapping parameter'
*********************************************************************
*       Initial values physical variables	
        w0	    'Abandoned derelict satellites 2023'		   /3500/
        z0	    'Rocket bodies 2023'				           /2050/
        f0	    'Number of fragments larger than 10cm 2023'	   /30000/
        d0	    'Total number of debris larger than 10cm 2023'	
        ss0	    'Operational satellites 2023'			       /8500/
        xx0	    'Satellites destroyed by collisions 2023'
        ff0	    'Total number of fragments larger than 1cm 2023'
        fff0 	'Total number of fragments larger than 1mm 2023'
        dd0	    'Total number of debris larger than 1cm 2023'
        ddd0 	'Total number of debris larger than 1mm 2023'
        ob1	    'Number of objects larger than 10cm (ESA)'	/30000/
        ob2	    'Number of objects larger than 1cm (ESA)'	/1000000/
        ob3	    'Number of objects larger than 1mm (ESA)'	/130000000/
        nu1	    'Ratio >1cm to >10cm'
        nu2	    'Ratio >1mm to >1cm'
        nu3	    'Ratio >1mm to >10cm'
	
*********************************************************************
**      Iteration Parameters
*********************************************************************
        max_iter    'Number of iteration'   /100/
        tolerance   'Tolerance'             /1e-8/
        ;
*********************************************************************


**	Calculus of key constants and initial values
    	mu	=	ss0/s0;
    	nu1 =  	ob2/ob1;
    	nu2 =  	ob3/ob2;
    	nu3 =  	ob3/ob1;
    	d0 	= 	w0+z0+f0;
        ff0	=	(1+nu1)*f0;
        dd0	= 	w0+z0+ff0;
        fff0=	nu3*f0;
        ddd0=	w0+z0+fff0;
        a0  =	y0/(k0**alpha1*s0**alpha2*N0**(1-alpha1-alpha2));
        x0	=	theta*dd0*s0;
        xx0 =  	theta*dd0*ss0;
*********************************************************************


*********************************************************************
*   	Timing
*********************************************************************
Parameters
    	beta(t)	'Discount factor'
    	tval(t)	'numerical value of t';
	    tval(t)	=	ord(t)-1;
*********************************************************************



*********************************************************************
*       TFP and ISTC trends
Parameters
*   Technology
        A(t)    'TFP'
        Q(t)    'Satellite ISTC'
        gA(t)   'Growth rate of TFG'
        gQ(t)   'Growth rate of Satellite ISTC'
        ;
    
        gA(t)           =   gA0*exp(-deltaA*(tval(t)-1));
        gQ(t)           =   gQ0*exp(-deltaQ*(tval(t)-1));
        a("2023")       =   a0;
        q("2023")       =   q0;
        loop(t,A(t+1)   =   a(t)*exp(gA(t)););
        loop(t,Q(t+1)   =   q(t)*exp(gQ(t)););
*********************************************************************



*********************************************************************
*       Launch cost trend
Parameters
        b(t)    'Fraction of launch cost'
        gb(t)   'Decline in the launch cost';
        
        gb(t)           =   gb0*exp(-deltab*(tval(t)-1));
        b("2023")       =   b0;
        loop(t,b(t+1)   =   b(t)*exp(gb(t)););
*********************************************************************        



*********************************************************************
*       Population
*********************************************************************
Parameters
        N(t)    'Population';
        N("2023")       =   N0;
        loop(t,N(t+1)   =   N(t)*(Nstar/N(t))**gN0;);
*********************************************************************
	


*********************************************************************
*       The terminal weight beta(tlast) computation
*********************************************************************
	   beta(tnotlast(t))	=	power(1+rho,-tval(t));
	   beta(tlast(t))		=	(1/rho)*power(1+rho,1-tval(t)-1);
*********************************************************************




*********************************************************************	
**      Exogenous block of space environment
*********************************************************************
Parameters
        s_trend(t)      'guess for satellites'
        ss_trend(t)     'guess for the bumber of satellites';
        s_trend(t)  =   s0*exp(0.01*(ord(t)-1));
        ss_trend(t) =   mu*s_trend(t);
*********************************************************************
*       Debris
*********************************************************************
Parameters
        ss(t)
        ll(t)
        h_trend(t)
        w(t)    'Guess of derelict satellites'
        z(t)    'Guess of rocket bodies'
        f(t)    'Guess of fragments larger than 10cm'
        ff(t)   'Guess of fragments larger than 1cm'
        d(t)    'Guess of debris larger than 10cm'
        dd(t)   'Guess of debris larger than 1cm'
        ;
        h_trend(t)=0.01;
        w(t)    =   w0*exp(0.01*(ord(t)-1));
        z(t)    =   z0*exp(0.01*(ord(t)-1));
        f(t)    =   f0*exp(0.01*(ord(t)-1));
        d(t)    =   w(t)+z(t)+f(t);
        ff(t)   =   (1+nu1)*f(t);
        dd(t)   =   w(t)+z(t)+ff(t);
        ss(t)   =   ss_trend(t);
        ll(t)   =   mu*h_trend(t)/eta;   
*********************************************************************
*       Damages
*********************************************************************        
Parameters
        x(t)    'Guess of value of satellite destroyed'
        xx(t)   'Guess of the number of satellite destroyed'
        ;
        x(t)    =   (1-v)*theta*dd(t)*s_trend(t);
        xx(t)   =   (1-v)*theta*dd(t)*ss_trend(t);
*********************************************************************



******************************************************************************************************
*           Variables of the decentralized economy
******************************************************************************************************
Variables
*           Economic variables of decentralized economy
	SU      'Total summatory discounted utility'
	y(t)	'Output'
   	c(t)    'Consumption'
   	rk(t)
   	rs(t)
   	wage(t)
   	ratcy(t)'Ratio consumption-output'
   	ik(t)	'Investment in Earth capital'
   	is(t)	'Investment in Space capital'
   	h(t)    'Investment in satellite'
   	l(t)    'Launch cost'
   	k(t)    'Earch capital stock'
	s(t)	'Stock of satellites'
	ii(t)   'Investment'
	hh(t)	'Number of new satellites inserted into orbit'
	U(t)	'Instantaneous utility'
	DU(t)   'Discounted utility'
	srk(t)	'Investment rate in Earth capital'
	srs(t)	'Investment rate in Space capital'
*    ll(t)   'Number of launches'
*ss(t)   'Number of operational satellites'
    ;
******************************************************************************************************
Equations
* 	         Economic variables
	utility(t)     'Instantaneous utility'
	sumdutility	   'Summatory of discounted utility'
	production(t)  'Cobb-Douglas production function'
	allocation(t)  'Household choose between consumption and investments'
	totalinvest(t) 'Total investment'
	launch(t)      'Launch cost'
	satel(t)       'Satellite cost'
	accumulak(t)   'Capital accumulation'
	accumulas(t)   'Satellite accumulation'
	ifinal(t)      'Minimal capital investment in final period'
	hfinal(t)      'Minimal satellite investment in final period'  
	iifinal(t)     'Minimal total investment in final period'
	ratiocy(t)     'Consumption-output ratio'
	ratioiy(t)     'Investment rate in Earth capital'
	ratiohy(t)     'Investment rate in Space capital'
*	nl(t)          'Number of launches'
*	ns(t)          'Number of satellites'
	ratek(t)
	rates(t)
	wages(t)
	limit(t)
	;
*******************************************************************************************************


********************************************************************************************************
**      Equations of the model
********************************************************************************************************
*       Utility
	utility(t)..	U(t)		=e=	((c(t)/N(t))**(1-sigma)-1)/(1-sigma);
	sumdutility..	SU  	 	=e=	sum(t, beta(t)*U(t)*N(t));
********************************************************************************************************
*       Production, consumption and investment
********************************************************************************************************
	production(t)..    		y(t) 		=e= 	a(t)*(k(t)**alpha1)*(s(t)**alpha2)*(N(t)**(1-alpha1-alpha2));
	allocation(t)..    		c(t) 		=e= 	rk(t)*k(t)+rs(t)*s(t)+wage(t)*N(t)-ik(t)-is(t);
	ratek(t)..              rk(t)       =e=     alpha1*y(t)/k(t);
	rates(t)..              rs(t)       =e=     alpha2*y(t)/s(t);
	wages(t)..              wage(t)     =e=     (1-alpha1-alpha2)*y(t)/N(t);
	totalinvest(t)..		ii(t)		=e=	    ik(t)+is(t);
	launch(t)..             l(t)        =e=     b(t)*is(t);
	satel(t)..              h(t)        =e=     is(t)-l(t);
	accumulak(tnotlast(t))..  	k(t+1) 		=e=	(1-deltak)*k(t)+ik(t);
	accumulas(tnotlast(t))..	s(t+1)		=e=	(1-deltas)*s(t)+q(t)*(1-b(t))*is(t)-x(t);
	ifinal(tlast)..     		ik(tlast) 	=g= 	(0.02+deltak)*k(tlast);
	hfinal(tlast)..			is(tlast)	=g=	(0.02+deltas)*s(tlast)/q(tlast);
	iifinal(tlast)..		ii(tlast)	=e=	ik(tlast)+is(tlast);
	ratiocy(t)..			ratcy(t)	=e=	c(t)/y(t);
	ratioiy(t)..			srk(t)		=e=	ik(t)/y(t);
	ratiohy(t)..			srs(t)		=e=	is(t)/y(t);
	limit(t)..                  x(t)        =l= s(t);


********************************************************************************************************
*       Bounds
********************************************************************************************************
	k.lo(t) = 0.001;
	s.lo(t)	= 0.001;
	c.lo(t) = 0.001;
	y.lo(t) = 0.001;
	is.lo(t)= 0;
********************************************************************************************************	


********************************************************************************************************
*       Initial conditions
********************************************************************************************************
	k.fx(tfirst) 	= 	k0;
	s.fx(tfirst) 	= 	s0;
	y.fx(tfirst)    =   y0;
********************************************************************************************************	


********************************************************************************************************
*       Solution options
********************************************************************************************************
   	option iterlim = 99900;
   	option reslim  = 99999;
   	option solprint = on;
	option limrow = 0;
	option limcol = 0;
********************************************************************************************************


********************************************************************************************************
*           Equation of the Model
********************************************************************************************************
model DISED / all /;
********************************************************************************************************


********************************************************************************************************
*               Temporary Variables
********************************************************************************************************
Parameters
    dso_temp(t)
    dsd_temp(t)
    rb_temp(t)
    rbd_temp(t)
    efrag_temp(t)
    ex_temp(t)
    eds_temp(t)
    erb_temp(t)
    edsb_temp(t)
    erbb_temp(t)
    x_old(t)
    xx_old(t)
    ;
********************************************************************************************************


********************************************************************************************************    
*               Loop iterations to update orbital debris
********************************************************************************************************
Set
    iter    /1*100/;
Parameter
    diff;
    diff=0.1;
Scalar
    iter_count;
********************************************************************************************************

    
********************************************************************************************************
*           Timing for the internal loop
********************************************************************************************************
Alias (t, tp);
Set tnext(t, tp) "Mapping from t to next period tp";
tnext(t, tp)$(ord(tp) = ord(t) + 1) = yes;
********************************************************************************************************



********************************************************************************************************
*           Main loop. Maximization with exogenous orbital debris (no internalization)
********************************************************************************************************
Loop(iter$(ord(iter) <= max_iter and diff > tolerance),
    iter_count=ord(iter);
    solve DISED maximizing SU using nlp;

** Update derelict satellite
    dso_temp(t)     =   xi*deltas*mu/q(t)*s.l(t);
** Update derelict satellites destroyed by collisions
    dsd_temp(t)     =   theta*(dd(t)+(1-v)*mu/q(t)*s.l(t))*w(t);
** Update rocket bodies
    rb_temp(t)      =   psi*ll(t);
** Update rocket bodies destroyed by collisions
    rbd_temp(t)     =   theta*(dd(t)+(1-v)*mu/q(t)*s.l(t))*z(t);
** Update fragments emissions by launches
    efrag_temp(t)   =   omega*ll(t);
** Update fragments emissions by operational satellite collisions
    ex_temp(t)      =   gammas*(1-v)*theta*dd(t)*mu/q(t)*s.l(t);
** Update fragments emissions by derelict satellite collisions
    eds_temp(t)     =   gammaw*theta*dd(t)*w(t);
** Update fragments emissions by rocket bodies collisions
    erb_temp(t)     =   gammaz*theta*dd(t)*z(t);
** Update fragments emissions by derelict satellite breakups
    edsb_temp(t)    =   phiw*epsilonw*w(t);
** Update fragments emissions by rocket bodies breakups
    erbb_temp(t)    =   phiz*epsilonz*z(t);
** Update dereclit satellites, rocket bodies, and fragments after each iteration
    loop(tnext(t,tp),
    w(tp)=(1-deltaw-epsilonw)*w(t)-dsd_temp(t)+dso_temp(t);
    z(tp)=(1-deltaz-epsilonz)*z(t)-rbd_temp(t)+rb_temp(t);
    f(tp)=(1-deltaf)*f(t)+efrag_temp(t)+ex_temp(t)+eds_temp(t)+erb_temp(t)+edsb_temp(t)+erbb_temp(t);
    );
** Update debris larger than 1 cm
    ff(t)   =   (1+nu1)*f(t);
    dd(t)   =   w(t)+z(t)+ff(t);
    ss(t)   =   mu/q(t)*s.l(t);
    ll(t)   =   mu*(1-b(t))*is.l(t)/eta;
** Store previous satellites destroyed to calculate the difference
    x_old(t) = x(t);
    xx_old(t)   =   xx(t);
** Update the damage function based on new population of debris
    x(t)=(1-v)*theta*dd(t)*s.l(t);
** Update the damage as number of operational satellites destroyed
    xx(t)=(1-v)*theta*dd(t)*ss(t);
** Calculate the difference between the current and previous satellites destroyed
    diff = smax(t, abs(xx(t)-xx_old(t)));
    
break$(diff <= tolerance);
);
********************************************************************************************************

* Export to Excel using GDX utilities
execute_unload "decentralized.gdx"
execute '=gdx2xls decentralized.gdx';
